Link del repositorio de github aquí

#cuando te metes
load("OUTPUT/objetos_actuales.Rdata")

#antes de salir
save.image(file="objetos_actuales.Rdata")

1. Importación de paquetes necesarios

library(tidyjson)
library(rjson)
library(tidyverse)
library(dplyr)
library(ggplot2)
library(lubridate)
library(readr)
library(plotly)
library(DT)

2. Introducción

En nuestro seminario vamos a estudiar las altas hospitalarias por ictus, en los diferentes hospitales de Castilla y León durante cuatro años, y como estas, están relacionadas con distintos factores de riesgo.

Lo primero de todo, explicar que los factores de riesgo son características o circunstancias que atentan contra el equilibrio, contra la salud, y que causan enfermedades o muerte. [senado1999factores]

Durante el desarrollo de nuestro seminario vamos a trabajar con datos de dos factores de riesgo, los cuales son: la calidad del aire y valorar si los centros de desintoxicación los podemos considerar como factor de riesgo o como factor de prevención de sufrir un ictus.

La calidad del aire es quizas un factor de riesgo del ictus, quizás poco conocido, pero estudios confirman que “Niveles más elevados de NO el día de ingreso conlleva un aumento del riesgo de mortalidad” [suarez2024influencia]

3. Objetivos

Objetivos principales:

1.Estudiar la cantidad de altas hospitalarias por ictus que hay en cada provincia de Castilla y León y ver como varían dependiendo del sexo y de la edad.

2.Analizar si a peor calidad del aire hay más riesgo de sufrir un ictus.

3.Estudiar la relación entre centros de ayuda a la desintoxicación por provincia y el número de ictus.

4. Datos utilizados

4.1 Procedencia de datos

Los datos han sido recabados de Datos abiertos del gobierno de españa que hasta el momento contiene cerca de 90000 conjuntos de datos, distribuidos en un amplio número de temáticas.También supera las 550000 distribuciones, Además de 313 iniciativas de datos abiertos.

4.2 Descripción y Importación de datos

Descripción:

Contamos con 1 set de datos en tipo .csv. Este formato, (Comma Separated Values) son archivos de texto que en pricipio deberían estar con los caracteres separados por comas, aunque en algunos de nuestros datos se usa punto y coma en lugar de la coma.

Este es nuestro archivo .csv.:

-calidad_aire

También contamos con dos sets de datos .json (JavaScript Object Notation) es un formato para almacenar e intercambiar datos.

Estos son nuestros archivos .json.:

-altas_hospitalarias_ictus

-centros_servicios_saisde

Este último proviene del Instituto Nacional de Estadística (INE).

4.3 Importación de datos

4.3.1 Importamos los datos de las altas hospitalarias por ictus

altas_hospitalarias_ictus <- fromJSON(file="INPUT/data/altas_hospitalarias_ictus.json")
#altas_hospitalarias_ictus

4.3.2 Importamos los datos de la calidad del aire

calidad_aire <- read_csv("INPUT/data/calidad_aire.csv")
#calidad_aire

Esta tabla decidimos no mostrarla porque es tan grande que no deja renderizarlo.

#datatable(calidad_aire)

4.3.3 Importamos los datos de los centros de desintoxicación

centros_servicios_saisde <- fromJSON(file = "INPUT/data/centros_servicios_saisde.json")
#centros_servicios_saisde

Importamos la función que cuenta el número de ictus por provincia:

Está función se creo antes de conocer tuberías, la función group_by y summarise con las cuales podemos obetener el mismo resultado pero lo dejaremos a modo comparativo, para ver las ventajas y eficiencia de trabajar con estas estructuras.

source("INPUT/functions/ictus_por_provincia.R")

Aquí vemos el ejemplo trabajando con dichas metodologías

prov_e_ictus <- Provincias%>%
  group_by(Provincia)%>%
  summarise(numero_ictus=n())

5. Manipulación de datos y resolución de cuestiones:

5.1 Ictus, edad y sexo

Para resolver el primer objetivo vamos a:

  • Obtener las edades de todos los pacientes para ver si la edad es un factor determinante en el aumento de los ictus, para ello necesitamos conocer como está organizado nuestro json.
altas_hospitalarias_ictus %>%
  spread_all() %>%
  gather_object %>%
  json_types %>%
  count(name, type)
## # A tibble: 4 × 3
##   name             type       n
##   <chr>            <fct>  <int>
## 1 datasetid        string 20090
## 2 fields           object 20090
## 3 record_timestamp string 20090
## 4 recordid         string 20090

Vemos que en fields tenemos un objeto con los datos que queremos extraer, por lo tanto hacemos:

altas_hospitalarias_ictus %>%
  enter_object(fields) %>%
  gather_object %>%
  spread_all %>%
  json_types%>%
  count(name, type)
## # A tibble: 11 × 3
##    name                                     type       n
##    <chr>                                    <fct>  <int>
##  1 ambito_de_procedencia                    string 20040
##  2 area                                     string 20090
##  3 dia_de_la_semana_en_la_fecha_del_alta    string 20090
##  4 dia_de_la_semana_en_la_fecha_del_ingreso string 20090
##  5 edad                                     number 20090
##  6 fecha_de_alta                            string 20090
##  7 fecha_de_ingreso                         string 20090
##  8 hospital                                 string 20090
##  9 provincia                                string 20090
## 10 sexo                                     string 20090
## 11 zona_basica_de_salud_del_paciente        string 20038

Dejamos la tabla entera tabulada, para trabajar a partir de ahora con este objeto:

altas_hosp_ictus_tab <- altas_hospitalarias_ictus%>%
  spread_all()%>%
  enter_object("fields")%>%
  spread_all()
#altas_hosp_ictus_tab
#datatable(altas_hosp_ictus_tab)

Debido a que esta tabla también es muy grande tampoco la mostramos con el paquete: library(DT). Debido a que esta tabla también es muy grande tampoco la mostramos con el paquete: library(DT).Debido a que esta tabla también es muy grande tampoco la mostramos con el paquete: library(DT).Debido a que esta tabla también es muy grande tampoco la mostramos con el paquete: library(DT).

Obtener el número de casos por provincias para empezar a relacionarlo con los siguientes objetivos.

Provincias <- altas_hosp_ictus_tab%>%
  rename(Provincia=fields.provincia)%>%
  select(Provincia)%>%
  mutate(Provincia=case_when(
    Provincia=="Ávila"~"Avila",
    .default = Provincia
  ))

#Provincias

Obtenemos el número de ictus por provincia

prov_e_ictus <- Provincias%>%
  group_by(Provincia)%>%
  summarise(numero_ictus=n())
datatable(prov_e_ictus)

Calculamos el número de personas que han sufrido un ictus para cada edad, para ello obtenemos el atributo edad de cada una de las personas. Además, lo haremos de dos maneras:

edades_total <- altas_hosp_ictus_tab %>%
  rename(Edad = edad) %>%         
  select(Edad)
#edades_total
#datatable(edades_total)

Esta tabla hemos decido no mostrarla puesto que es muy grande.

Importamos función que cuenta el total para cada una de las edades:

5.1.1 Primera manera (orientada en programación) para obtener el numero de casos por edad

Se trata de una función programada desde un punto de vista de programación, más adelante mostraremos como hacerlo desde un enfoque en el manejo de datos como se requiere en la asignatura y por ende, la gran mayoría de código quedará implementado de la segunda manera.

source("INPUT/functions/total_por_edad.R")
resultados_edad <- total_por_edad(edades_total)
datatable(resultados_edad)

5.1.2 Segunda manera (orientada al manejo de datos) para obtener el numero de casos por edad

ictus_por_edad <- edades_total%>%
  group_by(Edad)%>%
  summarise(num_ictus=n())
datatable(ictus_por_edad)

Podemos observar que de esta forma obtenemos menos filas, que de la primera forma, esto se debe a que en esta segunda tabla nos quita las personas con las edades que no han sufrido un ictus.

Regresión polinómica:

regresion_edad_ictus <- lm(num_ictus~poly(Edad,3),data=ictus_por_edad)
summary(regresion_edad_ictus)
## 
## Call:
## lm(formula = num_ictus ~ poly(Edad, 3), data = ictus_por_edad)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -227.10  -81.20  -21.57   79.80  279.25 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      202.93      10.82  18.759  < 2e-16 ***
## poly(Edad, 3)1  1449.03     107.64  13.462  < 2e-16 ***
## poly(Edad, 3)2  -527.86     107.64  -4.904 3.87e-06 ***
## poly(Edad, 3)3 -1368.52     107.64 -12.714  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 107.6 on 95 degrees of freedom
## Multiple R-squared:  0.7943, Adjusted R-squared:  0.7879 
## F-statistic: 122.3 on 3 and 95 DF,  p-value: < 2.2e-16

Hemos decidido realizar una regresión polinómica ya que una lineal no se ajusta exactamente lo que queremos y nos daba un R^2 muy bajo, en cambio la polinómica nos da un R^2 de 0.7994 lo cual nos indica que no es un mal modelo. Esto sucede porque en la regresión lineal debe crecer linealmente hasta el infinito pero a ciertas edades es biologicamente imposible que aumente ya que los humanos tenemos una vida limitada. En los resultados del modelor polinomial vemos que la edad es un factor estadísticamente siginificativo.

Gráfico que explica el modelo:

ggplot(data=ictus_por_edad, aes(x=Edad,y=num_ictus))+
  geom_point()+
  geom_smooth()+
  labs(x="Edad (años)",y="Número de ictus",title="Regresión polinómica edad-ictus")

Se ve claramente como a edades bajas el número deictus es muy reducido pero a partir de la barrera de los 45-50 empieza a subir drásticamente hasta los 75-88 que empieza a disminuir ya que a estas edades comienza a disminuir la tasa de supervivencia ya que las personas fallecen por otras patologías o por causas naturales.

Ahora vamos a obtener el número de hombres y mujeres para ver si el sexo influye en la cantidad de ictus

Seleccionamos el sexo:

sexototal <- altas_hosp_ictus_tab %>%
  rename(Sexo = sexo) %>%         
  select(Sexo)
total_sexo <- sexototal$Sexo
#datatable(sexototal)

Basado en métodos de programación

Llamamos a la función que nos cuenta el número de hombres y de mujeres que han sufrido un ictus:

source("INPUT/functions/total_por_sexo.R")
resultados_sexo <- total_por_sexo(total_sexo)
datatable(resultados_sexo)

Basado en técnicas de manejo de datos

ictus_por_sexo<- sexototal%>%
  group_by(Sexo)%>%
  summarise(num_ictus=n())%>%
  arrange(desc(Sexo))
datatable(ictus_por_sexo)

Un total de 11124 hombres sufrieron un ictus y el número de mujeres fue de 8966.

Gráfico de barras:

Mostramos de manera visual con un gráfico de barras sencillo la cantidad de hombres y mujeres que han sufrido un ictus:

ictus_por_sexo%>%
  ggplot(mapping=aes(x=Sexo, y=num_ictus))+
  geom_bar(stat="identity", aes(fill=Sexo))+
  scale_fill_manual(values=c("Hombre"="blue", "Mujer"="deeppink"))+
  labs(x="Sexo",y="Numero de ictus",title="Cantidad de personas (por sexo) que sufrieron un ictus")

ANOVA:

Realizamos un ANOVA para ver si el sexo influye en el riesgo de sufrir un ictus.

Para ello necesitamos conocer la media de la edad y su desviación estandar.

media_desviacion <- altas_hosp_ictus_tab %>%
  group_by(sexo) %>%
  summarise(
    media = mean(edad, na.rm = TRUE),
    desviacion_estandar = sd(edad, na.rm = TRUE)
  )

datatable(media_desviacion)
anova_resultado <- aov(edades_total$Edad ~ sexototal$Sexo, data = altas_hosp_ictus_tab)

summary(anova_resultado)
##                   Df  Sum Sq Mean Sq F value Pr(>F)    
## sexototal$Sexo     1  131314  131314   765.3 <2e-16 ***
## Residuals      20088 3446613     172                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Análisis estadístico:

Tomamos como Test Hipótesis:

H0: Medias iguales –> p-valor>0.05

H1: Medias distintas –>p-valor<0.05

Tras realizar el anova(analisis de la varianza) con la edad como variable numérica y el sexo como variable categórica vemos que el p-valor es muy bajo rechazamos la H0, lo que nos indica que el sexo es un factor estadisticamente significativo para los ictus.

5.2 Calidad del aire

Calculamos la media de los factores más destacados del aire, agrupados por fecha y provincia.

media_factores_aire <-calidad_aire%>%
    rename(provincia=Provincia)%>%
    mutate(Fechas = my(Fecha))%>%
    mutate(mes= month(Fechas))%>%
    mutate(anyo= year(Fechas))%>%
    group_by(Fechas, provincia,mes,anyo)%>%
    summarise(
      NO=mean(`NO (ug/m3)`, na.rm = TRUE),
      O3=mean(`O3 (ug/m3)`, na.rm = TRUE),
      PM25=mean(`PM25 (ug/m3)`, na.rm = TRUE)
      
    )%>%
    arrange(Fechas)
datatable(media_factores_aire)

Necesitamos cambiar el formato de la fecha, para así poder trabajar y agrupar datos que contienen fechas, ya que el paquete lubridate permite hacer lo que necesitamos con las fechas, en formato fecha.

Además, tenemos que usar un case_when para que Ávila pase a llamarse Avila para luego poder hacer diferentes operaciones necesarias.

altas_hosp_ictus_fecha_bien <- altas_hosp_ictus_tab %>%
  mutate(Fecha = ymd_hms(fields.fecha_de_ingreso))%>%
  mutate(mes= month(fields.fecha_de_ingreso))%>%
  mutate(anyo = year(fields.fecha_de_ingreso))%>%
  mutate(provincia=case_when(
    provincia=="Ávila"~"Avila",
    .default = provincia
  ))
#datatable(altas_hosp_ictus_fecha_bien)

Esta tabla no la mostramos tampoco debido a su gran tamaño.

Juntamos ambas tablas con sus atributos comunes que nos son de interes (mes, año y provincia)

aire_con_ictus <- full_join(altas_hosp_ictus_fecha_bien, media_factores_aire, by = c("mes","anyo", "provincia"))
  
aire_con_ictus_resumen <- aire_con_ictus%>%
  rename(ambito=fields.ambito_de_procedencia,dia_semana=fields.dia_de_la_semana_en_la_fecha_del_ingreso)%>%
  select(Fechas,mes,anyo,edad,sexo,dia_semana,provincia,NO,O3,PM25,ambito)
#datatable(aire_con_ictus_resumen)

Esta tabla tampoco la mostramos debido a su gran tamaño.

Calculamos la media de cada una de las sustancias (NO, O3 y PM25) agrupadas por provincia, año y mes. Además, contamos el número de ictus que hay para cada uno de ellos.

calcula_media_sustancias_por_provincia_y_meses <-aire_con_ictus_resumen %>%
  group_by(provincia, anyo, mes) %>%  
  summarize(num_ictus=n(),
            avg_NO = mean(NO, na.rm = TRUE),
            avg_O3=mean(O3, na.rm = TRUE),
            avg_PM25=mean(PM25, na.rm = TRUE)) 
datatable(calcula_media_sustancias_por_provincia_y_meses)

En esta tabla aparece representado la media de cada una de las sustancias del aire en cada provincia por mes y año.

Estudios confirman que “Niveles más elevados de NO el día de ingreso conlleva un aumento del riesgo de mortalidad debido a ictus” El gráfico que se muestra a continuación nos sirve para ver sí cuando la concentración de NO es mayor hay un mayor número de ictus.

calcula_media_sustancias_por_provincia_y_meses%>%
  ggplot(aes(x=avg_NO, y=num_ictus))+
  geom_point(aes(colour = factor(provincia)))+
  geom_smooth(method = "lm", aes(colour = factor(provincia)))+
  labs(x="Media concentración NO (ug/m3)", y="Número de ictus", title="Gráfico de dispersión del NO", colour="Provincias")

Podemos observar que las lineas para cada provincia según el articulo(ref) deberían salir con pendiente positiva pero no es así. Creemos que este resultado se debe a un pequeño matiz entre el artículo y nuestros datos. Dicho matiz se trata de que en el artículo nos indica que concentraciones altas de NO aumentaban el riesgo de mortalidad mientras que en nuestros datos los pacientes que tenemos han recibido el alta. Así que como no tenemos número de fallecidos no lo podemos contrastar realmente con el artículo.

El ozono es uno de los principales contaminantes del aire. Tal y como afirma el profesor Shaowei Wu, de la universidad Xi´an Jiaotong (China): “estudios previos han sugerido que la contaminación por ozono daña el corazón y los vasos sanguíneos”.

Por ello vamos a realizar un gráfico de dispersión en el que veamos como varía el número de ictus dependiendo de como sea la media de O3 en cada provincia.

Gráfico de dispersión del 03:

calcula_media_sustancias_por_provincia_y_meses%>%
  ggplot(aes(x=avg_O3, y=num_ictus))+
  geom_point(aes(colour = factor(provincia)))+
  geom_smooth(method = "lm", aes(colour = factor(provincia)))+
  labs(x="Media concentración O3 (ug/m3)", y="Número de ictus", title="Gráfico de dispersión del O3", colour="Provincias")

En este caso, sí que podemos observar una pendiente positiva a mayor concentración de ozono por lo que podemos concluir que si es un factor que aumenta el numero de ictus.

Las hospitalizaciones aumentan en un día determinado entre un 0,5% y un 1% por cada aumento de 10 mg/m3 de PM2.5. PM2.5: Hace referencia al material particulado respirable presente en la atmósfera. Se puede encontrar más información sobre PM2.5 aqui

Por ello vamos a realizar un gráfico de dispersión para ver si cuando los niveles de PM2.5 son más elevados hay mayor número de ictus.

Gráfico de dispersión del PM25:

calcula_media_sustancias_por_provincia_y_meses%>%
  ggplot(aes(x=avg_PM25, y=num_ictus))+
  geom_point(aes(colour = factor(provincia)))+
  geom_smooth(method = "lm", aes(colour = factor(provincia)))+
  labs(x="Media concentración PM2,5 (ug/m3)", y="Número de ictus", title="Gráfico de dispersión del PM2,5", colour="Provincias")

En este caso tenemos provincias que no tienen valores de PM2.5 (las que no tienen linea en el gráfico) pero en las que lo tenemos en todas se cumple que a mayor concentración mayor número de ictus.

Gráfico de barras de todos los gases:

Vamos a graficar todos los gases en una misma imagen, para ello, necesitamos modificar nuesta tabla aire_con_ictus_resumen añadiendo una columna llamada gases que contenga el nombre de los gases que hemos seleccionado y una columna media de concentración con los valores de dicha concentración:

tabla_con_columna_con_todos_gases <- aire_con_ictus_resumen%>%
  select(mes,anyo,provincia,NO,O3,PM25)%>%
  pivot_longer(names_to="gases",values_to="media_de_concentracion",cols=c(NO,O3,PM25))
#datatable(tabla_con_columna_con_todos_gases)

Faceta:

Realizamos una faceta para mostrar lo explicado anteriormente:

tabla_con_columna_con_todos_gases %>%
  ggplot(aes(x=anyo, y=media_de_concentracion))+
  geom_point(aes(colour = factor(provincia)))+
  geom_smooth(method = "lm", aes(colour = factor(provincia)))+
  facet_wrap(facets=vars(gases), nrow=1)+ 
  labs(x="Año",y="Media de concentración", title="Faceta de las sustancias del aire", colour="Provincias")

Ahora vamos a realizar 3 gráficos, uno para cada sustancia

Gráfico de barras de la media de NO:

Gráfico que muestra la media de NO (en una escala de colores) por provincia y mes:

calcula_media_sustancias_por_provincia_y_meses %>%
  ggplot(aes(x = reorder(factor(provincia),factor(anyo)), y = num_ictus)) + 
  geom_bar(stat = "identity", aes(fill = avg_NO)) +
  scale_fill_gradient(low = "turquoise", high = "blue" , space = "Lab", na.value = "grey50", guide = "colourbar") +
  labs(x = "", y = "Número ictus", fill = "Media de concentración de NO") +
  theme_classic()

En esta gráfica aparecen representados los números de ictus por provincia y la concentración de NO, la cual aparece pintada más oscura cuando es mayor y más clara cuando es menor.

p <- calcula_media_sustancias_por_provincia_y_meses %>%
  ggplot(aes(x = reorder(factor(provincia),factor(anyo)), y = num_ictus)) + 
  geom_bar(stat = "identity",aes(fill = avg_NO)) +
  scale_fill_gradient(low = "turquoise", high = "blue", space = "Lab", na.value = "grey50", guide = "colourbar") +
  labs(x = "Provincia y Año", y = "Número de ictus", fill = "Media de concentración de NO") +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))
p

# Convierte el gráfico ggplot a plotly
ggplotly(p, tooltip = c("x", "y", "fill"))

Esta gráfica es una versión mejorada de la anterior, en la cual aparecen representados los ictus que ha habido en cada provincia y en cada año, cada barra tiene 12 partes diferenciadas (una por mes) en las de los años 2020 en adelante; en el año 2019 aparecen menos ya que solo contamos con datos de octubre en adelante. Las barras están coloreadas de distintas tonalidades, corriespondiendose las oscuras a mayores concentraciones de NO y las claras a menores concentraciones de NO. Además pasando el cursor por las barras, gracias a la tecnología del paquete plotly, podemos ver podemos ver con que año y provincia se corresponde esa barra, la media de concentración de NO en ese mes y también el número total de ictus que han ocurrido ese mes.

Gráfico de barras de la media de O3:

Gráfico que muestra la media de 03 (en una escala de colores) por provincia y mes:

calcula_media_sustancias_por_provincia_y_meses %>%
  distinct(provincia, .keep_all = TRUE) %>%  
  ggplot(aes(x = reorder(provincia,anyo), y = num_ictus)) + 
  geom_bar(stat = "identity", aes(fill = avg_O3)) +
  scale_fill_gradient(low = "turquoise", high = "blue", space = "Lab", na.value = "grey50", guide = "colourbar") +
  labs(x = "", y = "Número ictus", fill = "Media de concentración de 03") +
  theme_classic()

En esta gráfica aparecen representados los números de ictus por provincia y la concentración de O3, la cual aparece pintada más oscura cuando es mayor y más clara cuando es menor.

p1 <- calcula_media_sustancias_por_provincia_y_meses %>%
  ggplot(aes(x = interaction(provincia, anyo), y = num_ictus, fill = avg_O3)) + 
  geom_bar(stat = "identity") +
  scale_fill_gradient(low = "turquoise", high = "blue", space = "Lab", na.value = "grey50", guide = "colourbar") +
  labs(x = "Provincia y Año", y = "Número de ictus", fill = "Media de concentración de O3") +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

# Convierte el gráfico ggplot a plotly
ggplotly(p1, tooltip = c("x", "y", "fill"))

Esta gráfica es una versión mejorada de la anterior, en la cual aparecen representados los ictus que ha habido en cada provincia y en cada año, cada barra tiene 12 partes diferenciadas (una por mes) en las de los años 2020 en adelante; en el año 2019 aparecen menos ya que solo contamos con datos de octubre en adelante. Las barras están coloreadas de distintas tonalidades, corriespondiendose las oscuras a mayores concentraciones de O3y las claras a menores concentraciones de O3. Además pasando el cursor por las barras, gracias a la tecnología del paquete plotly, podemos ver podemos ver con que año y provincia se corresponde esa barra, la media de concentración de NO en ese mes y también el número total de ictus que han ocurrido ese mes.

Gráfico de barras de la media de PM2,5:

Gráfico que muestra la media de PM2.5 (en una escala de colores) por provincia y mes:

calcula_media_sustancias_por_provincia_y_meses %>%
  distinct(provincia, .keep_all = TRUE) %>%  
  ggplot(aes(x = reorder(provincia,anyo), y = num_ictus)) + 
  geom_bar(stat = "identity", aes(fill = avg_PM25)) +
  scale_fill_gradient(low = "turquoise", high = "blue", space = "Lab", na.value = "grey50", guide = "colourbar") +
  labs(x = "", y = "Número ictus", fill = "Media de concentración de PM2.5") +
  theme_classic()

En esta gráfica aparecen representados los números de ictus por provincia y la concentración de O3, la cual aparece pintada más oscura cuando es mayor y más clara cuando es menor.

p2 <- calcula_media_sustancias_por_provincia_y_meses %>%
  ggplot(aes(x = interaction(provincia, anyo), y = num_ictus, fill = avg_PM25)) + 
  geom_bar(stat = "identity") +
  scale_fill_gradient(low = "turquoise", high = "blue", space = "Lab", na.value = "grey50", guide = "colourbar") +
  labs(x = "Provincia y Año", y = "Número de ictus", fill = "Media de concentración de PM2,5") +
  theme_classic() +
  theme(axis.text.x = element_text(angle = 90, hjust = 1))

# Convierte el gráfico ggplot a plotly
ggplotly(p2, tooltip = c("x", "y", "fill"))

Esta gráfica es una versión mejorada de la anterior, en la cual aparecen representados los ictus que ha habido en cada provincia y en cada año, cada barra tiene 12 partes diferenciadas (una por mes) en las de los años 2020 en adelante; en el año 2019 aparecen menos ya que solo contamos con datos de octubre en adelante. Las barras están coloreadas de distintas tonalidades, corriespondiendose las oscuras a mayores concentraciones de O3y las claras a menores concentraciones de O3. Además pasando el cursor por las barras, gracias a la tecnología del paquete plotly, podemos ver podemos ver con que año y provincia se corresponde esa barra, la media de concentración de NO en ese mes y también el número total de ictus que han ocurrido ese mes. Cabe destacar que en este caso en concreto, hay muchos valores NaN, que son los que aparecen de color gris.

Gráfico de barras de los días de ingreso:

Como dato curioso, queremos ver si el día de la semana está relacionado con los ictus

numero_ictus_dia_semana <-aire_con_ictus_resumen %>%
  group_by(dia_semana) %>%  
  summarize(num_ictus=n())
            
numero_ictus_dia_semana %>%
  ggplot(mapping=aes(x = dia_semana, y = num_ictus)) + 
  geom_bar(stat = "identity", fill="orchid") +
  labs(x = "Días de la semana", y="Número de ictus") +
  theme_classic()

Por lo leído y descrito en uno de los artículos, niveles altos de NO el día de ingreso conlleva a mas riesgo de mortalidad, al ser nuestros datos de altas hospitalarias por ictus y no de muertes, lo que vamos a intentar ver es si los días que hay un mayor número de ingresos coinciden con niveles altos de NO.

5.3 Centros de desintoxicación

Importamos los datos de los centros de desintoxicación.

Visionamos los datos y vemos que con hacer un spread_all podemos obtener lo que necesitamos para nuestro trabajo.

spread_all(centros_servicios_saisde)
Centros <- centros_servicios_saisde %>%
  spread_all() %>%
  mutate(fields.provincia=case_when(
    fields.provincia=="Ávila"~"Avila",
    .default = fields.provincia))%>%
  count(provincia = fields.provincia) %>%
  as.data.frame()
Centros
##    provincia  n
## 1      Avila  6
## 2     Burgos 12
## 3       León 18
## 4   Palencia 13
## 5  Salamanca 16
## 6    Segovia  5
## 7      Soria  4
## 8 Valladolid 17
## 9     Zamora  8
#datatable(Centros)

Importamos función que imprime los habitantes de castilla y león por provincias:

source("INPUT/functions/poblacion_cyl.R")
datatable(poblacion_cyl)

Para tener una referencia entre los centros de desintoxicación por provincia con el número de habitantes de esa provincia calcularemos cuál es el número de habitantes por centro; con el obetivo de ver realmente si hay más ictus en aquellas provincias donde el número de habitantes por centro es menor y el consumo de drogas es un factor significativo en los ictus o por el contrario no lo es.

pob_centros <- full_join(poblacion_cyl,Centros)
tabla <- pob_centros%>%
  select(provincia,Poblacion_total,n)
source("INPUT/functions/divide.R")
hab_centro <- divide(tabla$Poblacion_total,tabla$n)
#faltaria indicar que con que provincia se corresponden esos datos.
tabla$hab_centro <- hab_centro
resultado <- tabla %>%
  select(provincia, hab_centro)
datatable(resultado)

Lo lógico es esperar que en las provincias que cuentan con un menor número de habitantes, haya menor número de ictus y menor número de centros de desintoxicación.

numero_ictus <- ictus_por_provincia(Provincias)
df <- as.data.frame(numero_ictus)



pob_prov_y_total <- poblacion_cyl%>%
  select(provincia, Poblacion_total)

datatable(pob_prov_y_total )
pob_cyl_sin_ult_fila <- pob_prov_y_total [1:9,]
datatable(pob_cyl_sin_ult_fila)
arrange(.data=Centros,(Centros$provincia))
##    provincia  n
## 1      Avila  6
## 2     Burgos 12
## 3       León 18
## 4   Palencia 13
## 5  Salamanca 16
## 6    Segovia  5
## 7      Soria  4
## 8 Valladolid 17
## 9     Zamora  8
datos <- data.frame(pob=pob_cyl_sin_ult_fila,
                    centros= arrange(.data=Centros,(Centros$provincia)),
                    num_ictus=df$numero_ictus)

datatable(datos)

Gráfico de barras en el que aparece representado el número de centros por provincia:

datos%>%
   ggplot(mapping = aes(x = reorder(pob.provincia, num_ictus), y = centros.n)) +
  geom_bar(stat = "identity",fill="blue") +
  labs(x = "Provincias", y = "Número de centros") +
  theme_classic()

Gráfica de barras en la que se representa el número de ictus en cada provincia:

Cabe destacar que se colorean de un color más fuerte las barras con más centros de desintoxicación.

datos %>%
  ggplot(mapping = aes(x = reorder(pob.provincia, num_ictus), y = num_ictus)) +
  geom_bar(stat = "identity", aes(fill = centros.n)) +
  scale_fill_gradient(low = "peachpuff", high = "red") +
  labs(x = "Provincias", y = "Número de ictus", fill = "Centros detox") +
  theme_classic() 

Gráfico de barras en el que representamos el número de habitantes por centro para cada provincia:

resultado
##          provincia hab_centro
## 1            Avila   26627.33
## 2           Burgos   29780.83
## 3             León   24920.72
## 4         Palencia   12137.46
## 5        Salamanca   20443.06
## 6          Segovia   31066.40
## 7            Soria   22382.00
## 8       Valladolid   30666.65
## 9           Zamora   20865.88
## 10 Castilla y León         NA
resultado_sin_ult_fila=resultado[1:9,]




ggplot(mapping = aes(x = reorder(datos$pob.provincia, datos$num_ictus), y = resultado_sin_ult_fila$hab_centro)) +
  geom_bar(stat = "identity", fill = "green") +
  labs(x = "Provincias", y = "Número habitantes por centro") +
  theme_classic() 

Obtenemos una tabla que nos permite visualizar de manera numérica toda la información estudiada anteriormente:

tabla_centros <- tibble(Provincia=datos$pob.provincia,                     
                            Numero_de_ictus=datos$num_ictus,                                          Poblacion_total_cyl=pob_cyl_sin_ult_fila$Poblacion_total,
                            Numero_centros=datos$centros.n,
                        Habitantes_por_centro=resultado_sin_ult_fila$hab_centro,
                        Porcentaje_ictus_por_hab=(Numero_de_ictus/Poblacion_total_cyl)*100)
datatable(tabla_centros)
# Calcular métricas adicionales
tabla_metricas <- tabla_centros %>%
  mutate(tasa_incidencia = (Numero_de_ictus / Poblacion_total_cyl) * 1000,
         densidad_centros = (Numero_centros / Poblacion_total_cyl) * 10000)
ggplot(tabla_metricas, aes(x = Numero_centros, y = Numero_de_ictus)) +
  geom_point() +
  geom_smooth() +
  labs(title = "Relación entre numero de centros y numero de ictus",
       x = "Centros s",
       y = "Número de ictus ")

modelo_regresion_polinomico <- lm(
  Numero_de_ictus ~ poly(Numero_centros,4),
  data = tabla_metricas
)
summary(modelo_regresion_polinomico)
## 
## Call:
## lm(formula = Numero_de_ictus ~ poly(Numero_centros, 4), data = tabla_metricas)
## 
## Residuals:
##        1        2        3        4        5        6        7        8 
##  -229.86   977.06  -333.56 -1180.64   -10.78   -11.33    89.28   586.70 
##        9 
##   113.14 
## 
## Coefficients:
##                          Estimate Std. Error t value Pr(>|t|)   
## (Intercept)                2232.2      282.7   7.895  0.00139 **
## poly(Numero_centros, 4)1   3547.4      848.2   4.182  0.01389 * 
## poly(Numero_centros, 4)2    633.9      848.2   0.747  0.49637   
## poly(Numero_centros, 4)3    183.0      848.2   0.216  0.83972   
## poly(Numero_centros, 4)4   -839.9      848.2  -0.990  0.37815   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 848.2 on 4 degrees of freedom
## Multiple R-squared:  0.8267, Adjusted R-squared:  0.6533 
## F-statistic: 4.769 on 4 and 4 DF,  p-value: 0.07972

https://fhernanb.github.io/libro_regresion/predict.html

tabla_predicciones <- data.frame(tabla_metricas, predicciones = predict(modelo_regresion_polinomico))

# Visualización del impacto
ggplot(tabla_predicciones, aes(x = Habitantes_por_centro, y = tasa_incidencia)) +
  geom_point() +
  geom_line(aes(y = predicciones), color = "blue") +
  labs(title = "Impacto estimado de centros detox en la tasa de ictus",
       x = "Habitantes por centro",
       y = "Tasa de Ictus por 1,000 habitantes")

Tras la realización de las tres gráficas anteriores para ver como están relacionados el número de ictus de cada provincia con el número de centros de desintoxicación que hay en las mismas, y el número de habitantes por centro de cada una de ellas, hemos llegado a las siguientes conclusiones:

  • Las poblaciones que cuentan con un menor número de habitantes (Ávila, Palencia, Soria y Zamora) son a su vez las que cuentan con un menor porcentaje de ictus; esto seguramente pueda deberse a que por regla general en las ciudades pequeñas se suelen llevar mejores habitos de vida.

  • Podemos destacar tres casos en especial:

    1. En Palencia, cuentan con un número elevado de centros de desintoxicación para la población que tiene; lo cual supone que haya pocos habitantes por centro, pudiendo ver así como el porcentaje de ictus es el segundo más pequeño.Por lo tanto podría parecer que la presencia de más centros ayuda a reducir el número de ictus.

    2. Segovia tiene un gran número de habitantes por centro, el más alto de todos, lo cual podría indicar que las personas que acuden van a recibir un “peor” tratamiento ya que el tratamiento se va a tener que repartir entre más pacietes y al ser peor el riesgo de ictus va a ser mayor, sin embargo en esta población no se cumple la hipotésis planteada en la conclusión de Palencia.

    3. Salamanca es la provincia que cuenta con un mayor porcentaje de ictus por habitante, lo cual es destacable ya que no es la que más habitantes tiene y además es la segunda que menos habitantes por centro tiene.

Por lo tanto,visualmente y falta de un estudio estadístico preciso no podemos determinar si el número de centros es un factor estadísticamente significativo para el riesgo de padecer un ictus.

Conclusiones

Queda demostrado que tanto la edad como el sexo influyen en el riesgo de sufrir un ictus, esto lo hemos comprobado realizando una regresión polinómica (edad) y un ANOVA (sexo), con sus correspondientes gráficas para hacerlo más visual.

En cuanto a las sustancias del aire, podemos decir que la concentración de NO y la de PM2,5 si que se relacionan con la probabilidad de sufrir un ictus; mientras que la concentración de O3 no podemos llegar a saber si verdaderamente esta relacionada ya que lo que está comprobado es que el hecho de que haya elevadas concentraciones de NO aumenta la mortalidad por ictus en ese día, y nuestros datos no son de muertes sino de altas hospitalarias.

En relación con la cantidad de centros, tampoco podemos determinar si el que haya más numero de centros reduce la probabilidad de sufrir un ictus.

.